function [gu,nXu,ic,binary_waveInPool] = empirical_density_generate_gu_cubic(Xr_wave,Xu_pooled,gr_wave)
% This function generates the g for unique households, and
% count the number of unique households in two-dimensional space. 
% This is used for making consistent covariate space across waves, because 
% not all entries in Xu_pooled are present in Xr_wave

% Inputs: 
% (1) Xr_wave: sorted covariates in wave
% (2) Xu_pooled: unique X in the pooled sample
% (3) gr_wave: sorted g in wave

% Outputs: 
% (1) gu: sum of gr at the unique household level
% (2) nXu: number of households at each unique household level
% (3) ic: index mapping to unique household for each household in wave
% (4) binary_waveInPool: check whether wave unique household is a subset of
%                        unique households in the pooled sample

%% 

[Xu_wave,ia,ic] = unique(Xr_wave,'rows','stable'); 
gu_wave = arrayfun(@(x) sum(gr_wave(ic==x)),1:size(Xu_wave,1),'uni',1)';
nXu_wave=accumarray(ic,1); 

% check where each row (iR) of "Xu_wav"e is in "Xu_pooled", 
% sum this row (iR) in "gu_wave" into corresponding row in "gu_pooled"
gu=zeros(size(Xu_pooled,1),1); nXu = gu;
binary_waveInPool = ismember(Xu_pooled,Xu_wave,'rows');
gu(binary_waveInPool)=gu_wave;
nXu(binary_waveInPool)=nXu_wave;

end

